Introduction

In this document, we will outline the Bayesian analogs of the statistical analyses described here (and in the previous course lecture).

Helper functions. This function extracts results from different models and generate results of the same format to be used in visualizations

tidy.wrapper = function(model) {
  if (class(model) == "lm") {
    tidy(model, conf.int = TRUE) %>%
      select(-c(statistic, p.value)) %>%
      mutate(model = "Frequentist") %>%
      select(model, everything())
  } else if (class(model) == "brmsfit") {
    tidy(model) %>%
      filter(effect == "fixed") %>%
      select(c(term:conf.high)) %>%
      mutate(model = "Bayesian") %>%
      select(model, everything())
  } else {
    stop("unknown model class")
  }
}

Dataset 1

Let’s first load and take a look at the data in a table:

dataset = readr::read_csv("data/blinded.csv")

head(dataset)
experiment condition effectiveness
1 no_graph 7
1 graph 6
1 no_graph 9
1 no_graph 4
1 graph 6
1 no_graph 7

Recall from the previous lecture that this dataset contains results from four experiments. In each experiment, participants are placed in either a graph or no graph condition, and asked to rate the effectiveness of the intervention on a nine-item Likert scale. Here, we plot the data to show the distribution of these responses:

dataset %>% 
  mutate(effectiveness = fct_rev(factor(effectiveness, levels = 1:9)),
         experiment = as.factor(experiment)) %>%
  
  # stacked bar plot
  ggplot(aes(x = condition, fill = effectiveness)) +
  geom_bar(position = "stack", stat="count") +
  
  # plot data for different experiments as small multiples
  facet_wrap(. ~ experiment) +
  
  # grey color scale is robust for colorblind
  scale_fill_brewer(palette="Purples", drop = FALSE) +
  
  # horizontal plot
  coord_flip() +
  
  # legend
  guides(fill = guide_legend(reverse = TRUE))

For the purposes of this lecture, we will confine ourselves to the first experiment.

exp1.data = dataset %>%
  filter(experiment == 1)

head(exp1.data)
experiment condition effectiveness
1 no_graph 7
1 graph 6
1 no_graph 9
1 no_graph 4
1 graph 6
1 no_graph 7

Model 1. Wilcoxon signed rank test

The first model discussed in the previous lecture was the Wilcoxon signed rank test. This is a non-parametric test which we will skip for now. Although, there exists Bayesian non-parametric methods, they are more advanced for this lecture.

Model 2. Students’s t test

The second model discussed was the T-test. Although R contains a function (t.test) to perform this analysis, the t-test is essentially a linear regression, and thus can be performed using the linear model function in R (lm). This is the linear model equivalent for the paired sample t-test:

Frequentist t-test

dataset1.lm.freqt <-
  lm(
    effectiveness ~ condition - 1, # we remove intercept
    data = exp1.data
  )

Step 1: Build the model (likelihood)

Before we implement a Bayesian model, let us take a look at the distribution of responses. Although we know that we are going to use a t distribution, we still plot the data to set how it looks like. Usually, this will help us decide our likelihood function.

  exp1.data %>%
     ggplot(., aes(x = effectiveness)) +
  geom_histogram(fill = DARK_PURPLE, color = NA, binwidth = 0.5, center = 0) +
  scale_x_continuous(breaks = seq(0, 10, by = 1))

Let’s also look how a student t distribution varies different parameter. Mu is location, sigma is dipsersion, and nu (df) control the tail.

ggplot() +  xlim(-10, 10) +
  geom_function(
    color = "black",
    fun = function(x)
      dstudent_t(x, mu = 0, sigma = 1, df = 5) ,
    size = 1
  ) +
  geom_function(
    color = "gray80",
    fun = function(x)
      dstudent_t(x, mu = 0, sigma = 1, df = 50),
    size = 1
  ) +
  geom_function(
    color = "gray60",
    fun = function(x)
      dstudent_t(x, mu = 0, sigma = 4, df = 5),
    size = 1
  ) +
  geom_function(
    color = "gray40",
    fun = function(x)
      dstudent_t(x, mu = 2,  sigma = 1, df = 5),
    size = 1
  ) +
  xlab("x")

We can see from the above plot that the responses are discrete. This is important to keep in mind when specifying a Bayesian model. Next, we will implement the Bayesian analog of the linear model described earlier. Here’s what the model formula will look like when implemented using the brm formula (bf function):

dataset1.brm.bayesiant.formula <- bf(effectiveness ~ condition - 1,
                                     family = student())

Let us take a minute to understand this code. It says that we are regressing the variable condition on effectiveness. The - 1 removes the intercept term. The family argument is used to specify the probability distribution of the likelihood — in other words, what is the distribution of \(P(y)\).

Step 1: Build the model (priors)

The blank ones are flat (uniform) priors. These are improper priors and usually needed to be set. We can and should adjust other priors given by brm.

as_tibble(get_prior(dataset1.brm.bayesiant.formula, data = exp1.data))
prior class coef group resp dpar nlpar bound source
b default
b conditiongraph default
b conditionno_graph default
gamma(2, 0.1) nu default
student_t(3, 0, 2.5) sigma default

Next, we need to specify priors for the parameters in the model. We will discuss briefly how these prior distributions were obtained.

For the prior on b which is the mean parameter of the student t distribution (used as the likelihood), a unbiased assumption would be that they should be centered around 5 (which is the center of the 9-item likert scale), and the mean would likely be greater than 1 and less than 9 (unless every single participant responded either 1 or 9; but we know that was not the case).

The prior on sigma determines the standard deviation parameter of the student t distribution. Now, considering that our data is bounded between 1 and 9, a uniform distribution will have the maximum variance given these constraints. The variance of such an uniform distribution would be: sd(runif(1e4, 1, 9) which is approximately 2.5. We know that our data will most likely have less variance than such a uniform distribution. Thus, we want the prior on sigma to assign very little probability mass for values greater than 2.5. Our prior student_t(3, 0, 1) assigns less than 0.05 probability to values greater than 2.5 (you can check by running the following code in the console: sum(gamlss.dist::rTF(1e4, 0, 1, 3) > 2.5) / 1e4).

dataset1.brm.bayesiant.priors = c(
      prior(normal(5, 1.5), class = "b"),
      prior(student_t(3, 0, 1), class = "sigma")
  )

Step 1: Build the model (prior predictive checks)

Before we implement the regression model, it is advisable to perform some prior predictive checks. Bayesian models are generative — in other words, it is possible to sample values from the prior distributions, feed them through the likelihood function to obtain a prior predictive distribution.

The prior predictive distribution should ideally resemble the data generating process, and should not assign probability mass to unlikely or impossible values. If the prior predictive distribution is assigning substantial probability mass to unlikely or impossible values, we should want to adjust our choice of prior distributions. Prior predictive checks also help make explicit some of the assumptions that go into the prior specification process.

Another important thing to note is that brms often specifies improper priors (denoted by flat) by default. However, we cannot sample draws from the prior predictive distribution if the priors are improper.

The following code block implements some prior predictive checks:

dataset1.brm.bayesiant.priorchecks <-
  brm(
    dataset1.brm.bayesiant.formula,
    prior = dataset1.brm.bayesiant.priors,
    data = exp1.data,
    backend = BRM_BACKEND,
    sample_prior = "only",
    file = "rds/dataset1.brm.bayesiant.priorchecks.rds"
  )

n_prior_draws <- 30

# extract n = 100 draws from the prior predictive distribution
dataset1.bayesiant.yprior <-
  posterior_predict(dataset1.brm.bayesiant.priorchecks, ndraws = n_prior_draws)

# the following computes the probability density for each "draw"
dataset1.bayesiant.ppc = dataset1.bayesiant.yprior %>%
  as_tibble() %>%
  mutate(.draw = row_number()) %>%
  pivot_longer(
    cols = starts_with("V"),
    names_to = "participant",
    names_prefix = "V",
    values_to = ".value"
  ) %>%
  group_by(.draw) %>%
  summarise(.value = list(.value)) %>%
  mutate(dens.x = map(.value, ~ density(.)$x),
         dens.y = map(.value, ~ density(.)$y))

dataset1.bayesiant.ppc %>%
  select(-.value) %>%
  unnest(c(dens.x, dens.y)) %>%
  ggplot() +
  geom_line(
    aes(x = dens.x, y = dens.y, group = .draw),
    color = DARK_PURPLE,
    alpha = .5,
    size = .75
  ) +
  labs(y = "", x = 'Draws from the prior predictive distribution')

# Fumeng -> Abhraneel: I rewrote it...
# I actually don't immediately know how to implement density estimation per .draw

# flat it
dim(dataset1.bayesiant.yprior) <- n_prior_draws * 123

# create a new table and assign .draw number
tibble(.value = dataset1.bayesiant.yprior,
       .draw = rep(1:n_prior_draws, times = 123)) %>%
  ggplot() +
  geom_density(
    aes(x = .value, group = .draw),
    color = DARK_PURPLE,
    alpha = .5,
    size = .5
  )

Although, based on the above plot, it seems like the model is assigning some non-zero probability to values which are not possible (y <= 0 and y > 9), these are primarily a result of our choice of the likelihood function — the student t distribution does not allow us to set bounds on the predictive distribution. Another reason for extreme values may be the use of a student’s t prior for sigma with 3 degrees of freedom — this distribution has fat tails and thus might result in predicting large values for the standard deviation parameter in our model.

For now, let’s move forward, and fit the model.

Step 2: Computing the posterior probability distribution

Next, we compute the posterior probability distribution using brms and Stan. Depending on the complexity of your model, this step may take a lot of time. Our model is not so complex, hence this will not take too long :)

dataset1.brm.bayesiant =
  brm(
    dataset1.brm.bayesiant.formula,
    prior = dataset1.brm.bayesiant.priors,
    data = exp1.data,
    backend = BRM_BACKEND,
    # There are many other parameters you can play with.
    # Here we use default parameters to simplify the lecture.
    # We save the model 
    # But if you have this file, brm only loads the model from this file
    file = "rds/dataset1.brm.bayesiant.rds"
  )

Step 3: Evaluate the fit

After the compilation step of the model is complete, we can evaluate the fit and examine the result. Similar to lm(), the first step is to call the summary() function on the variable which contains the model compilation result:

summary(dataset1.brm.bayesiant)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: effectiveness ~ condition - 1 
##    Data: exp1.data (Number of observations: 123) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## conditiongraph        6.61      0.18     6.25     6.96 1.00     3840     2692
## conditionno_graph     6.77      0.17     6.43     7.11 1.00     3211     2355
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.30      0.11     1.07     1.52 1.00     3172     2490
## nu       16.84     11.62     4.36    46.74 1.00     3081     2832
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

This will give us an overview of the coefficients of the various parameters, along with the 95% credible intervals. This result shows that there is substantial overlap between the credible intervals for the two conditions. Before we delve deeper into the results, it is important to run some diagnostics on the sampling process.

Step 3: Evaluate posterior predictions

First, we call bayesplot to draw posterior distributions and MCMC traces for each parameter. We want to check whether the four chairs have mixed well. If they have, this implies that the model has converged.

color_scheme_set(scheme = "purple")
plot(dataset1.brm.bayesiant, newpage = T)

Next, we need to examine if the posterior distributions computed by our MCMC process resembles the data. Each draw from the posterior predictive distribution should roughly resemble a “hypothetical experiment” with the same experimental design. Thus, if our posterior predictive distribution looks very different from the data, it implies that something has gone awry in our model building step (Step 1).

They are many ways you can draw from the posterior predictive distribution. Here we use posterior_predict from brms.

dataset1.bayesiant.y <- exp1.data$effectiveness

dataset1.bayesiant.yrep <-
  posterior_predict(dataset1.brm.bayesiant)

dataset1.bayesiant.yrepgroup <-
  posterior_predict(dataset1.brm.bayesiant, data.frame(condition = c('graph', 'no_graph')))

head(as_tibble(dataset1.bayesiant.yrep[,1:11]))
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
6.440270 5.977234 6.900874 4.701375 7.170452 3.263761 6.592275 8.581487 4.573106 5.016169 7.093037
6.673812 8.328171 5.944403 7.329556 7.573959 4.559968 7.088508 7.700899 7.527388 5.979694 9.111188
6.680615 6.707617 9.736384 9.343441 8.087803 4.691732 5.436256 9.025326 5.613129 6.243677 5.763645
7.165171 3.938078 8.055621 7.618785 6.061229 4.580845 6.698171 3.310805 4.854779 4.683616 7.132660
5.600135 4.455032 6.197594 5.334783 8.090987 7.699026 6.578442 7.225055 6.790332 4.850451 6.231267
5.400997 5.004001 4.820987 7.282429 7.165793 7.385508 5.583039 8.057038 6.371572 5.565271 7.011921

We use bayesplot to perform these diagnostic tests. First, we plot the first eight draws of posterior predictions and the original data. Again, there are many ways to do this (and these will be discussed in greater detail in Lecture 3).

ppc_hist(y = dataset1.bayesiant.y,
         yrep = dataset1.bayesiant.yrep[1:8,],
         binwidth = .5)

An alternative approach would be to plot densities for each draw from the posterior predictive distribution superimposed with the density from the actual data:

ppc_dens_overlay(y = dataset1.bayesiant.y,
                 yrep = dataset1.bayesiant.yrep[1:30,])

Finally, we look at whether our model is able to capture the summary statistics properly. Below we show the distribution of the estimated posterior mean for the two conditions, as well as the actual mean of the data in the two conditions. Since the actual mean appears to be centered around the estimated posterior distribution of the mean, we can conclude that it is roughly capturing the relevant information from the experimental data.

ppc_stat_grouped(y = dataset1.bayesiant.y,
                 yrep = dataset1.bayesiant.yrep,
                 group = exp1.data$condition, binwidth = .1)

Step 3: Comparing posterior distributions

Now that we have performed model diagnostics, we would like to interpret the results. The pertinent research question here is whether there is a difference in the effectiveness rating between the graph and no graph condition. We use conditional distributions to compare the mean difference between the two conditions. add_epred_draws() from tidybayes is the most convenient way.

dataset1.bayesiant.posterior_epred <-
  tibble(condition = c('graph', 'no_graph')) %>%
  add_epred_draws(dataset1.brm.bayesiant,
                  re_formula = NA,
                  allow_new_levels = FALSE) %>%
  ungroup()

head(dataset1.bayesiant.posterior_epred)
condition .row .chain .iteration .draw .epred
graph 1 NA NA 1 6.475753
graph 1 NA NA 2 6.745635
graph 1 NA NA 3 6.589682
graph 1 NA NA 4 6.552249
graph 1 NA NA 5 6.525029
graph 1 NA NA 6 6.525103

Here add_epred_draws() takes two arguments:

  • newdata: is used to generate predictions. This variable should describe the experimental design. In this example, our experimental design consists of two conditions, graph and no graph, which is being passed as the first argument.

  • object: the model fit object

We then transform the data, subtract the two conditions from each other (using the compare_levels() function), and calculate the credible intervals (using mean_qi()).

dataset1.bayesiant.posterior_comparison <-
  dataset1.bayesiant.posterior_epred %>%
  select(-c(.chain, .iteration, .row)) %>%
  compare_levels(variable = .epred, by = condition)

head(dataset1.bayesiant.posterior_comparison)
.draw condition .epred
1 no_graph - graph 0.1617545
2 no_graph - graph 0.1336147
3 no_graph - graph 0.2767903
4 no_graph - graph 0.2008134
5 no_graph - graph 0.0610346
6 no_graph - graph 0.0427875 
dataset1.bayesiant.posterior_comparison %>%
        mean_qi(.epred)
condition .epred .lower .upper .width .point .interval
no_graph - graph 0.1579494 -0.3061581 0.6378437 0.95 mean qi

A basic way to interpret and communicate the results is to plot them. We plot this result using credible interval. Lecture 3 will cover more ways of representing uncertainty for effective communication.

dataset1.bayesiant.posterior_epred %>%
  ggplot(aes(x = .epred, fill = condition)) +
  #geom_point(aes(x = .epred, group = condition), color = DARK_PURPLE, size = 3) +
  geom_density(
    #fill = DARK_PURPLE,
    color = NA,
    alpha = .5,
    #geom = 'dots',
    #size = 1,
    #height = 0
  ) +
  scale_fill_manual(values = COLOR_PALETTE) + 
  #geom_vline(aes(xintercept = 0), linetype = "dashed", color = "black") +
  #coord_cartesian(xlim = c(-1, 1))  +
  xlab('')

dataset1.bayesiant.posterior_comparison %>%
  ggplot(aes(x = .epred)) +
  #geom_point(aes(x = .epred, group = condition), color = DARK_PURPLE, size = 3) +
  geom_density(
    fill = DARK_PURPLE,
    color = NA,
    alpha = .5,
    #geom = 'dots',
    #size = 1,
    #height = 0
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "black") +
  coord_cartesian(xlim = c(-1, 1))  +
  ggtitle('no_graph - graph') + 
  xlab('difference in mean')

dataset1.bayesiant.posterior_comparison %>%
  median_qi(.epred) %>%
  ggplot() +
  geom_point(aes(x = .epred, y = condition), color = DARK_PURPLE, size = 3) +
  geom_errorbarh(
    aes(xmin = .lower, xmax = .upper, y = condition),
    color = DARK_PURPLE,
    alpha = .5,
    size = 2,
    height = 0
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(0, 2), xlim = c(-1, 1))  +
  xlab('') + ylab('')

dataset1.bayesiant.posterior_comparison %>%
  mutate(if_greater = .epred < 0) %>%
  summarise(`Pr(graph > no_graph)` = mean(if_greater))
condition Pr(graph > no_graph)
no_graph - graph 0.25975

Finally, we compare the results obtained from the Bayesian model to the frequentist estimates:

bind_rows(tidy.wrapper(dataset1.lm.freqt),
          tidy.wrapper(dataset1.brm.bayesiant)) %>%
  ggplot() +
  geom_pointrange(
    aes(
      color = model,
      y = estimate,
      ymin = conf.low,
      ymax = conf.high,
      x = term
    ),
    position = position_dodge(width = 0.2)
  ) +
  scale_color_brewer(palette = "Set1") +
  ylab('effectiveness') +
  scale_y_continuous(breaks = 1:9, limits = c(1, 9)) +
  coord_flip()

Model 3. Ordinal Logistic Regression

A more appropriate to analyze ordianl data is to use ordinal logistic regression. In this section, we reanalyze the data using Bayesian ordinal logistic regression.

Step 1: Build the model (likelihood)

As always the first step is to specify the appropriate model formula:

dataset1.brm.olr.formula <-
  bf(effectiveness ~ condition,
     family = cumulative("logit"))

Before we go further into the model building step, it is important to understand what a ordinal regression model is estimating. The figure below shows the cumulative proportion of participants who responded at least \(k\), i.e. \(Pr(y_i \leq k)\), on the likert questionnaire.

# figure from: Statistical rethinking with brms, ggplot2, and the tidyverse: Second edition by A Solomon Kurz
p1 = exp1.data %>%
  count(effectiveness) %>%
  mutate(pr_k = n / nrow(exp1.data),
         cum_pr_k = cumsum(pr_k)) %>%
  ggplot(aes(x = effectiveness, y = cum_pr_k,
             fill = effectiveness)) +
  geom_line(color = "#b8925f") +
  geom_point(colour = "#b8925f",
             size = 2.5,
             stroke = 1) +
  scale_x_continuous("Effectiveness", breaks = 1:9) +
  scale_y_continuous("Cumulative Proportion",
                     breaks = seq(0, 1, by = 0.2),
                     limits = c(0, 1)) +
  theme(axis.ticks = element_blank(),
        legend.position = "none")

p2 = exp1.data %>%
  count(effectiveness) %>%
  mutate(pr_k = n / nrow(exp1.data),
         log_cum_odds = logit(cumsum(pr_k))) %>%
  ggplot(aes(x = effectiveness, y = log_cum_odds,
             fill = effectiveness)) +
  geom_line(color = "#b8925f") +
  geom_point(colour = "#b8925f",
             size = 2.5,
             stroke = 1) +
  scale_x_continuous("Effectiveness", breaks = 1:9) +
  scale_y_continuous(
    "Log-Cumulative Odds",
    breaks = seq(-5, 5, by = 1),
    limits = c(logit(0.005), logit(0.99))
  ) +
  theme(axis.ticks = element_blank(),
        legend.position = "none")

ggarrange(p1, p2)

Since cumulative proportion is a value bounded between 0 and 1, linear regression models apply the logit function to apply the following transformation \(f: [0, 1] \to (-Inf, Inf)\). Applying this transformation to the data gives the plot on the right.

Ordinal regression estimates this cumulative probability: \(Pr(y_i \leq k)\). They do this using a series of parameters \(\alpha_k\) where \(k \in \{1, 2, 3, 4, 5, 6, 7, 8, 9\}\), according to the following equation:

\[ \begin{align} log(\frac{Pr(y_i <= k)}{1 - Pr(y_i <= k)}) = \alpha_k \end{align} \] Ordinal regression estimates this cumulative probability: \(Pr(y_i \leq k)\). Once we have estimated the cumulative probability, we can take the difference between successive items, \(Pr(y_i = k) = Pr(y_i \leq k) - Pr(y_i \leq k - 1)\) to estimate the discrete probability for each outcome:

# primary data
exp1.data_plot = exp1.data %>%
  count(effectiveness) %>%
  mutate(pr_k = n / nrow(exp1.data)) %>%
  add_row(effectiveness = 2, n = 0, pr_k = 0) %>%
  arrange(effectiveness) %>%
  mutate(cum_pr_k = cumsum(n / nrow(exp1.data))) %>% 
  mutate(discrete_probability = ifelse(effectiveness == 1, cum_pr_k, cum_pr_k - pr_k))

text = exp1.data_plot %>%
  mutate(
    text = effectiveness,
    effectiveness = effectiveness + 0.25,
    cum_pr_k = ifelse(cum_pr_k - 0.05 < 0.065, 0.05, cum_pr_k - 0.05)
  )

exp1.data_plot %>% 
  ggplot(aes(x = effectiveness, y = cum_pr_k)) +
  geom_line(aes(color = cum_pr_k), color = "#b8925f") +
  geom_linerange(aes(ymin = 0, ymax = cum_pr_k), alpha = 1/2, color = "#b8925f") +
  geom_linerange(aes(x = effectiveness,
                     ymin = ifelse(effectiveness == 1, 0, discrete_probability), 
                     ymax = cum_pr_k),
                 color = "black") +
  geom_point(colour = "#b8925f", size = 2.5, stroke = 1) +
  # number annotation
  geom_text(data = text, 
            aes(label = text),
            size = 4) +
  scale_x_continuous("Effectiveness", breaks = 1:9) +
  scale_y_continuous("Cumulative Proportion", breaks = seq(0, 1, by = 0.2), limits = c(0, 1)) +
  theme(axis.ticks = element_blank(),
        axis.title.y = element_text(angle = 90),
        legend.position = "none")

Step 1: Build the model (priors)

The next step is to specify the prior distributions. Let us first look at the default prior distributions:

as_tibble(get_prior(dataset1.brm.olr.formula, data = exp1.data))
prior class coef group resp dpar nlpar bound source
b default
b conditionno_graph default
student_t(3, 0, 2.5) Intercept default
Intercept 1 default
Intercept 2 default
Intercept 3 default
Intercept 4 default
Intercept 5 default
Intercept 6 default
Intercept 7 default
Intercept 8 default

We notice that the difference parameter, b does not have a proper prior, so we should specify one. The Intercept term does have a prior, but since the distributions are going to be logit transformed it is difficult to intuitively determine whether these are strong or diffuse prior distributions. This is another important question that can be answered through prior predictive checks:

# priors
dataset1.brm.olr.priors = c(
          prior(normal(0, 1), class = "b")
)

dataset1.brm.olr.priorchecks <- brm(
    dataset1.brm.olr.formula,
    data = exp1.data,
    prior = dataset1.brm.olr.priors,
    iter = 15000,
    warmup = 7500,
    backend = BRM_BACKEND,
    sample_prior = 'only',
    file = "rds/dataset1.brm.bayesiant.priorchecks1.rds"
  )

dataset1.olr.yprior <-  posterior_predict(dataset1.brm.olr.priorchecks, ndraws = 100)

ggplot() +
  geom_histogram(aes(x = dataset1.olr.yprior),
               fill = DARK_PURPLE,
               alpha = .5, size = 1,
               binwidth = 0.5, center = 0) +
 scale_x_continuous(breaks = 1:9, limits = c(0.5, 9.5)) + 
 labs(x = 'Prior predictive distribution',  y = "") +
 theme(
   axis.text.y = element_blank()
 )

Our prior predictive check reveals that the model is assigning most of the probability mass to middle items in the Likert scale. This is not necessarily a bad thing, and it might reflect how participants behave. However, it is possible that these priors are still diffuse, as noticed by the low number of Bulk_ESS and Tail_ESS. These may act as another diagnostic measure, and ideally you’d want these values to be at least greater than 1000. Let’s us try a different prior

dataset1.brm.olr.priors = c(prior(normal(0, 1), class = "b"),
                            prior(student_t(3, 0, 1.5), class = "Intercept"))

dataset1.brm.olr.priorchecks <- brm(
    dataset1.brm.olr.formula,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    iter = 15000,
    warmup = 7500,
    backend = BRM_BACKEND,
    sample_prior = 'only',
    file = "rds/dataset1.brm.bayesiant.priorchecks2.rds"
    )

summary(dataset1.brm.olr.priorchecks)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
##   Draws: 4 chains, each with iter = 15000; warmup = 7500; thin = 1;
##          total post-warmup draws = 30000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -3.23      2.22    -9.39    -0.25 1.01      757      182
## Intercept[2]         -1.68      1.45    -4.88     0.42 1.00      947      360
## Intercept[3]         -0.83      0.97    -2.82     0.97 1.00    19222    18896
## Intercept[4]         -0.26      0.91    -2.07     1.50 1.00    10572    18752
## Intercept[5]          0.26      0.91    -1.49     2.09 1.00     3883    10659
## Intercept[6]          0.83      0.99    -0.97     2.91 1.00     3430     1194
## Intercept[7]          1.62      1.24    -0.45     4.49 1.00     6057     1876
## Intercept[8]          3.40      3.13     0.21    10.12 1.00     5841     2686
## conditionno_graph    -0.01      1.00    -1.96     1.95 1.00    13088    12542
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Changing the priors appear to improve Bulk_ESS and Tail_ESS slightly. Let’s take a look at the prior predictive distribution:

dataset1.olr.yprior <-
  posterior_predict(dataset1.brm.olr.priorchecks, ndraws = 100)


ggplot() +
  geom_histogram(aes(x = dataset1.olr.yprior),
               fill = '#351c75',
               alpha = .5,
               size = 1,
               binwidth = .5) +
  scale_x_continuous(breaks = 1:9, limits = c(0.5, 9.5)) + 
  xlab('prior draws') + ylab('')

Although these priors may be improved, it is often a case of trial and error. For now, let us proceed with these priors.

Step 2: Computing the posterior probability

After we’ve completed the model building phase, we can then compile the model.

dataset1.brm.olr1 =
  brm(
    dataset1.brm.olr.formula,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    backend = BRM_BACKEND,
    file = "rds/dataset1.brm.olr1.rds"
  )

summary(dataset1.brm.olr1)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.79      0.94    -6.96    -3.23 1.00     1762     1141
## Intercept[2]         -4.20      0.78    -5.98    -2.93 1.00     2466     1865
## Intercept[3]         -3.46      0.57    -4.71    -2.47 1.00     3340     2834
## Intercept[4]         -2.24      0.34    -2.97    -1.60 1.00     4238     3087
## Intercept[5]         -1.32      0.26    -1.83    -0.81 1.00     4691     3695
## Intercept[6]         -0.30      0.23    -0.74     0.17 1.00     4264     3276
## Intercept[7]          1.28      0.27     0.77     1.81 1.00     3594     3170
## Intercept[8]          2.26      0.34     1.63     2.96 1.00     3702     2956
## conditionno_graph     0.19      0.30    -0.38     0.78 1.00     4122     2809
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Notice that we get a few error messages about divergent transitions. Here we adjust a few sampling parameters, namely adapt_delta and max_treedepth which are passed to the control argument, to help the MCMC sampler and avoid divergent transitions:

dataset1.brm.olr2 =
  brm(
    dataset1.brm.olr.formula,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    backend = BRM_BACKEND,
    warmup = 1500,
    iter = 2500,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "rds/dataset1.brm.olr2.rds"
  )

summary(dataset1.brm.olr2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
##   Draws: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.86      1.08    -7.47    -3.21 1.00     2304     1698
## Intercept[2]         -4.23      0.83    -6.17    -2.90 1.00     2907     2007
## Intercept[3]         -3.47      0.57    -4.75    -2.47 1.00     3605     3118
## Intercept[4]         -2.24      0.35    -2.95    -1.59 1.00     3962     3378
## Intercept[5]         -1.32      0.27    -1.86    -0.81 1.00     4105     3240
## Intercept[6]         -0.30      0.24    -0.75     0.17 1.00     4242     3291
## Intercept[7]          1.28      0.26     0.79     1.79 1.00     4248     3274
## Intercept[8]          2.25      0.33     1.63     2.94 1.00     4093     3365
## conditionno_graph     0.19      0.30    -0.41     0.78 1.00     4069     2912
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Step 3: Evaluate the fit

We then perform the same model diagnostics steps that we had performed for the previous model.

plot(dataset1.brm.olr2, newpage = T)

Step 3: Evaluate the posterior predictions

dataset1.olr.y <- exp1.data$effectiveness
dataset1.olr.yrep <- posterior_predict(dataset1.brm.olr2)

ppc_hist(y = dataset1.olr.y,
         yrep = dataset1.olr.yrep[1000:1007, ],
         binwidth = .5)

ppc_dens_overlay(y = dataset1.olr.y,
                 yrep = dataset1.olr.yrep[2000:2050, ], adjust = .5)

ppc_bars(y = dataset1.olr.y,
        yrep = dataset1.olr.yrep, binwidth = .5)

ppc_stat_grouped(y = dataset1.olr.y,
                 yrep = dataset1.olr.yrep,
                 group = exp1.data$condition)

Step 3: Compare posterior distributions

Finally we estimate the difference between the two conditions and plot the results:

dataset1.olr.posterior_epred <-
  tibble(condition = c('graph', 'no_graph')) %>%
  add_epred_draws(dataset1.brm.olr2,
                  re_formula = NA,
                  allow_new_levels = FALSE) %>%
  ungroup()

head(dataset1.olr.posterior_epred)
condition .row .chain .iteration .draw .category .epred
graph 1 NA NA 1 1 0.0136202
graph 1 NA NA 2 1 0.0144437
graph 1 NA NA 3 1 0.0102796
graph 1 NA NA 4 1 0.0056930
graph 1 NA NA 5 1 0.0272409
graph 1 NA NA 6 1 0.0041625 
dataset1.olr.posterior_epred  %>% 
  filter(.draw == 100) %>%
  group_by(condition) %>%
  rename(`probability` = .epred) %>%
  ggplot(., aes(x = .category, y = `probability`, fill = condition, color = condition)) + 
  geom_bar(stat = 'identity', alpha = .5, color = NA, width = .3) + 
  geom_point(aes(group = condition)) + 
  geom_path(aes(group = condition)) + 
  facet_wrap(condition ~ ., ncol = 2) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .2)) + 
  scale_color_manual(values = c(DARK_PURPLE, GOLD_COLOR)) + 
  scale_fill_manual(values = c(DARK_PURPLE, GOLD_COLOR)) +
  xlab('')

dataset1.olr.posterior_epred %>% 
  select(-.chain, -.iteration, -.row) %>%
  group_by(condition, .draw, .category) %>%
  rename(probability = .epred) %>%
  ungroup() %>%
  group_by(.category, condition) %>%
  median_qi(probability) %>%
ggplot(., aes(x = .category, y = probability, fill = condition, color = condition)) + 
  geom_bar(stat = 'identity', alpha = .5, color = NA, width = .3) + 
  geom_errorbar(aes(ymin = .lower, ymax = .upper), width = 0, size = 1) + 
  geom_point(aes(group = condition)) + 
  geom_path(aes(group = condition)) + 
  facet_wrap(condition ~ ., ncol = 2) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .2)) + 
  scale_color_manual(values = c(DARK_PURPLE, GOLD_COLOR)) + 
  scale_fill_manual(values = c(DARK_PURPLE, GOLD_COLOR))  +
  xlab('')

dataset1.olr.posterior_epred  %>% 
  filter(.draw == 100) %>%
  #arrange(.category) %>%
  group_by(condition) %>%
  mutate(`cumulative probability` = cumsum(.epred)) %>%
  ggplot(., aes(x = .category, y = `cumulative probability`, fill = condition, color = condition)) + 
  geom_bar(stat = 'identity', alpha = .5, color = NA, width = .3) + 
  geom_point(aes(group = condition)) + 
  geom_path(aes(group = condition)) + 
  facet_wrap(condition ~ ., ncol = 2) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .2)) + 
  scale_color_manual(values = c(DARK_PURPLE, GOLD_COLOR)) + 
  scale_fill_manual(values = c(DARK_PURPLE, GOLD_COLOR))  +
  xlab('')

dataset1.olr.posterior_epred %>% 
  select(-.chain, -.iteration, -.row) %>%
  group_by(condition, .draw) %>%
  arrange(.category) %>%
  mutate(`cumulative probability` = cumsum(.epred)) %>%
  ungroup() %>%
  group_by(.category, condition) %>%
  median_qi(`cumulative probability`) %>%
ggplot(., aes(x = .category, y = `cumulative probability`, fill = condition, color = condition)) + 
  geom_bar(stat = 'identity', alpha = .5, color = NA, width = .3) + 
  geom_errorbar(aes(ymin = .lower, ymax = .upper), width = 0, size = 2) + 
  geom_point(aes(group = condition)) + 
  geom_path(aes(group = condition), size = 1) + 
  facet_wrap(condition ~ ., ncol = 2) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .2)) + 
  scale_color_manual(values = c(DARK_PURPLE, GOLD_COLOR)) + 
  scale_fill_manual(values = c(DARK_PURPLE, GOLD_COLOR))  +
  xlab('')

dataset1.olr.posterior_comparison <-
  dataset1.olr.posterior_epred %>%
  select(-c(.chain, .iteration, .row)) %>%
  group_by(.category) %>%
  compare_levels(variable = .epred, by = condition)

head(dataset1.olr.posterior_comparison %>%
       mean_qi())
.category condition .epred .lower .upper .width .point .interval
1 no_graph - graph -0.0019084 -0.0124921 0.0052925 0.95 mean qi
2 no_graph - graph -0.0011591 -0.0087963 0.0034439 0.95 mean qi
3 no_graph - graph -0.0024396 -0.0153921 0.0071908 0.95 mean qi
4 no_graph - graph -0.0098904 -0.0462694 0.0218328 0.95 mean qi
5 no_graph - graph -0.0137213 -0.0614789 0.0312567 0.95 mean qi
6 no_graph - graph -0.0150607 -0.0680144 0.0333193 0.95 mean qi 
dataset1.olr.posterior_comparison %>%
  mean_qi(.epred) %>%
  ggplot() +
  geom_point(aes(y = .epred, x = .category), size = 3, color = DARK_PURPLE) +
  geom_errorbar(
    aes(ymin = .lower, ymax = .upper, x = .category),
    width = 0,
    size = 2,
    color = DARK_PURPLE,
    alpha = .5
  ) +
  coord_cartesian(ylim = c(-.1, .1)) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
  xlab('') + ylab('no_graph - graph') +
  theme(axis.title.y = element_text(angle = 0, vjust = .5))

dataset1.olr.posterior_comparison  %>%
  ggplot() +
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
  geom_density(
    aes(y = .epred),
    color = NA, 
    fill = DARK_PURPLE,
    alpha = .5
  ) +
  coord_cartesian(ylim = c(-.1, .1)) + 
  xlab('') + ylab('no_graph - graph') +
  theme(axis.title.y = element_text(angle = 0, vjust = .5),
        axis.text.x  = element_blank())

dataset1.olr.posterior_comparison  %>%
       ungroup() %>%
       mean_qi(.epred) %>%
  ggplot() +
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
  geom_point(aes(y = .epred, x = 1), size = 3, color = DARK_PURPLE) +
  geom_errorbar(
    aes(ymin = .lower, ymax = .upper, x = 1),
    width = 0,
    size = 2,
    color = DARK_PURPLE,
    alpha = .5
  ) +
  coord_cartesian(xlim = c(0, 2), ylim = c(-.1, .1)) + 
  xlab('') + ylab('no_graph - graph') +
  theme(axis.title.y = element_text(angle = 0, vjust = .5),
        axis.text.x  = element_blank())

Dataset 2

Load the data.

dataset2 = readr::read_csv("data/exp2.csv") %>%
  mutate(condition = condition == 'expansive') %>%
  group_by(participant)

head(dataset2)
participant adjP10 adjP20 adjP30 adjPDiff condition gender index change subgroup orig BIS
1 56.50000 55.80000 58.40000 1.90000 TRUE female 69 3.362832 expansiveHigh 56.90000 high
2 27.87500 40.14286 36.00000 8.12500 FALSE female 60 29.147982 constrictiveLow 34.67262 low
3 61.00000 59.00000 76.50000 15.50000 TRUE female 55 25.409836 expansiveLow 65.50000 low
4 24.57143 32.75000 37.85714 13.28571 FALSE female 79 54.069767 constrictiveHigh 31.72619 high
5 64.71429 54.00000 41.00000 -23.71429 TRUE female 69 -36.644592 expansiveHigh 53.23810 high
6 35.14286 52.40000 45.60000 10.45714 FALSE female 43 29.756098 constrictiveLow 44.38095 low

We plot the data to give us a better picture about its distribution.

dataset2 %>%
  mutate(c = as.factor(condition)) %>%
  ggplot(aes(x = change)) +
  geom_histogram(
    aes(y = ..density..),
    binwidth = 10,
    fill = DARK_PURPLE,
    alpha = .5,
    color = 'white'
  ) +
  geom_density(size = 1,
               adjust = 1.5,
               color = '#9281bf') +
   geom_function(
    color = "#222222",
    linetype = 'dashed',
    fun = function(x)
      dstudent_t(x, mu = 16,  sigma = 39, df = 6),
    size = 1
  ) + 
  scale_x_continuous(limits = c(-200, 200))

Step 1: Build the model (likelihood and priors)

This model is the BEST test model as described by Kruschke in the paper Bayesian estimation supersedes the t-test. In this model, \(\beta\) indicates the mean difference in the outcome variable between the two groups (in this case, the percent change in the BART scores). We fit different priors on \(\beta\) and set different weights on these priors to obtain our posterior estimate.

\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \alpha_{0} + \beta * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta &\sim \mathrm{N}(\mu_{0}, \sigma_{0}) \\ \sigma_a, \sigma_b &\sim \mathrm{Cauchy}(0, 2) \\ \nu &\sim \mathrm{exp}(1/30) \end{align} \]

We translate the above specification to R code.

dataset2.brm.student.formula <- bf(change ~ condition,
                                   sigma ~ condition,
                                   family = student())

head(as_tibble(get_prior(dataset2.brm.student.formula, data = dataset2)))
prior class coef group resp dpar nlpar bound source
b default
b conditionTRUE default
student_t(3, 22.6, 38.8) Intercept default
gamma(2, 0.1) nu default
b sigma default
b conditionTRUE sigma default 
# This breaks :-)
dataset2.brm.student.priorchecks = brm(
  dataset2.brm.student.formula,
  prior = c(
    prior(normal(0, 2), class = "b"),
    prior(cauchy(0, 2), class = "b", dpar = "sigma"),
    prior(exponential(0.033), class = "nu"),
    prior(student_t(3, 0, 5), class = "Intercept"),
    prior(student_t(3, 0, 1), class = "Intercept", dpar = "sigma")
  ),
  data = dataset2,
  backend = BRM_BACKEND,
  sample_prior = 'only',
  file = "rds/dataset2.brm.student.priorchecks.rds"
)
min(dataset2$change)
## [1] -36.64459
max(dataset2$change)
## [1] 200.6135
dataset2.yprior <-  posterior_predict(dataset2.brm.student.priorchecks)

ggplot() +
  geom_histogram(aes(x = dataset2.yprior),
               fill = DARK_PURPLE,
               alpha = .5, 
               size = 1, 
               center = 0) +
 #scale_x_log10() + 
 # scale_x_continuous(breaks = 1:9, limits = c(0.5, 9.5)) + 
 labs(x = 'Prior predictive distribution',  y = "") +
 theme(
   axis.text.y = element_blank()
 )

dim(dataset2.yprior) <- 4000 * 80

quantile(dataset2.yprior, probs = c(.025, .975))
##      2.5%     97.5% 
## -565.4208  614.8581

Step 2: Computing the posterior probability

dataset2.brm.student = brm(
 bf(change ~ condition, 
 sigma ~ condition,
 family = student()),
 prior = c(
   prior(normal(0, 2), class = "b"),
   prior(cauchy(0, 2), class = "b", dpar = "sigma"),
   prior(exponential(.033), class = "nu"),
   prior(student_t(3, 0, 5), class = "Intercept"),
   prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
 ), 
 data = dataset2,
 backend = BRM_BACKEND,
 file = "rds/dataset2.brm.student.rds"
)

dataset2.brm.student$prior
prior class coef group resp dpar nlpar bound source
normal(0, 2) b user
b conditionTRUE default
cauchy(0, 2) b sigma user
b conditionTRUE sigma default
student_t(3, 0, 5) Intercept user
student_t(3, 0, 2) Intercept sigma user
exponential(0.033) nu user

Step 3: Evaluate the fit (summary)

summary(dataset2.brm.student)
##  Family: student 
##   Links: mu = identity; sigma = log; nu = identity 
## Formula: change ~ condition 
##          sigma ~ condition
##    Data: dataset2 (Number of observations: 80) 
##   Draws: 4 chains, each with iter = 1000; warmup = 0; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              20.95      4.80    12.02    30.97 1.00     2901     2206
## sigma_Intercept         3.36      0.20     2.98     3.75 1.00     2534     2195
## conditionTRUE          -0.19      1.92    -3.98     3.55 1.00     3894     2821
## sigma_conditionTRUE     0.14      0.22    -0.30     0.58 1.00     4263     2970
## 
## Family Specific Parameters: 
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu     4.56      5.29     1.61    14.40 1.01     1753     1433
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Step 3: Evaluate the posterior predictions

plot(dataset2.brm.student, newpage = T)

dataset2.student.y <- dataset2.brm.student$data$change
dataset2.student.yrep <- posterior_predict(dataset2.brm.student)

dataset2.student.epred <- tibble(condition = c(TRUE, FALSE)) %>%
  add_epred_draws(dataset2.brm.student,
                  re_formula = NA,
                  allow_new_levels = TRUE) %>%
  ungroup()
ppc_hist(y = dataset2.student.y,
         yrep = dataset2.student.yrep[100:107, ],
         binwidth = 10)

ppc_dens_overlay(y = dataset2.student.y,
                 yrep = dataset2.student.yrep[2000:2030, ])

ppc_stat_grouped(
  y = dataset2.student.y,
  yrep = dataset2.student.yrep,
  group = dataset2$condition,
  binwidth = 5
)

Step 3: Compare the posterior predictions

dataset2.student.epred_comparison <-
  tibble(condition = c(TRUE, FALSE)) %>%
  add_epred_draws(dataset2.brm.student,
                  re_formula = NA,
                  allow_new_levels = FALSE) %>%
  ungroup() %>%
  select(-c(.chain, .iteration, .row)) %>%
  compare_levels(variable = .epred, by = condition) %>%
  rename(diff = .epred)


head(dataset2.student.epred_comparison)
.draw condition diff
1 TRUE - FALSE 0.0054167
2 TRUE - FALSE 0.2985710
3 TRUE - FALSE 1.1666700
4 TRUE - FALSE 1.4400100
5 TRUE - FALSE 0.5096380
6 TRUE - FALSE -1.4863400 
dataset2.student.epred_comparison %>%
  select(diff) %>%
  mean_qi() %>%
  ggplot() +
  geom_point(aes(x = diff, y = condition), size = 3, color = DARK_PURPLE) +
  geom_errorbarh(
    aes(xmin = .lower, xmax = .upper, y = condition),
    height = 0,
    color = DARK_PURPLE,
    alpha = .5,
    size = 2
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(0, 2), xlim = c(-5, 5))  +
  scale_y_discrete(label = c('expansive - not expansive')) +
  xlab('')  + ylab('')

Compare variance

dataset2.student.sigma.epred_comparison <-
  tibble(condition = c(TRUE, FALSE)) %>%
  add_epred_draws(dataset2.brm.student,
                  dpar = "sigma",
                  re_formula = NA,
                  allow_new_levels = FALSE) %>%
  ungroup() %>%
  select(-c(.chain, .iteration, .row)) %>%
  compare_levels(variable = sigma, by = condition) %>%
  rename(sigma.diff = sigma)


head(dataset2.student.sigma.epred_comparison)
.draw condition sigma.diff
1 TRUE - FALSE 0.5262710
2 TRUE - FALSE 6.1196667
3 TRUE - FALSE 5.9907197
4 TRUE - FALSE 9.7824273
5 TRUE - FALSE -0.5522292
6 TRUE - FALSE -5.9265989 
dataset2.student.sigma.epred_comparison %>%
  select(sigma.diff) %>%
  mean_qi() %>%
  ggplot() +
  geom_point(aes(x = sigma.diff, y = condition), size = 3, color = DARK_PURPLE) +
  geom_errorbarh(
    aes(xmin = .lower, xmax = .upper, y = condition),
    height = 0,
    color = DARK_PURPLE,
    alpha = .5,
    size = 2
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(0, 2), xlim = c(-15, 15))  +
  scale_y_discrete(label = c('expansive - not expansive')) +
  xlab('')  + ylab('')

Packages Information

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cmdstanr_0.4.0.9001 knitr_1.37          bayesplot_1.9.0    
##  [4] ggdist_3.1.1        tidybayes_3.0.2     brms_2.16.3        
##  [7] Rcpp_1.0.8.2        modelr_0.1.8        broom.mixed_0.2.7  
## [10] broom_0.7.12        gtools_3.9.2        forcats_0.5.1      
## [13] tidyr_1.2.0         ggpubr_0.4.0        ggplot2_3.3.5      
## [16] purrr_0.3.4         tibble_3.1.6        dplyr_1.0.8        
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3     ggsignif_0.6.3       ellipsis_0.3.2      
##   [4] ggridges_0.5.3       markdown_1.1         base64enc_0.1-3     
##   [7] rstudioapi_0.13      farver_2.1.0         rstan_2.21.3        
##  [10] bit64_4.0.5          svUnit_1.0.6         DT_0.21             
##  [13] fansi_1.0.2          mvtnorm_1.1-3        diffobj_0.3.5       
##  [16] bridgesampling_1.1-2 codetools_0.2-18     splines_4.1.3       
##  [19] shinythemes_1.2.0    jsonlite_1.8.0       shiny_1.7.1         
##  [22] readr_2.1.2          compiler_4.1.3       backports_1.4.1     
##  [25] assertthat_0.2.1     Matrix_1.4-0         fastmap_1.1.0       
##  [28] cli_3.2.0            later_1.3.0          htmltools_0.5.2     
##  [31] prettyunits_1.1.1    tools_4.1.3          igraph_1.2.11       
##  [34] coda_0.19-4          gtable_0.3.0         glue_1.6.2          
##  [37] reshape2_1.4.4       posterior_1.2.1      carData_3.0-5       
##  [40] jquerylib_0.1.4      vctrs_0.3.8          nlme_3.1-155        
##  [43] crosstalk_1.2.0      tensorA_0.36.2       xfun_0.30           
##  [46] stringr_1.4.0        ps_1.6.0             mime_0.12           
##  [49] miniUI_0.1.1.1       lifecycle_1.0.1      rstatix_0.7.0       
##  [52] zoo_1.8-9            scales_1.1.1         vroom_1.5.7         
##  [55] colourpicker_1.1.1   hms_1.1.1            promises_1.2.0.1    
##  [58] Brobdingnag_1.2-7    parallel_4.1.3       inline_0.3.19       
##  [61] RColorBrewer_1.1-2   shinystan_2.6.0      yaml_2.3.5          
##  [64] gridExtra_2.3        loo_2.5.0            StanHeaders_2.21.0-7
##  [67] sass_0.4.0           stringi_1.7.6        highr_0.9           
##  [70] dygraphs_1.1.1.6     checkmate_2.0.0      pkgbuild_1.3.1      
##  [73] rlang_1.0.2          pkgconfig_2.0.3      matrixStats_0.61.0  
##  [76] distributional_0.3.0 evaluate_0.15        lattice_0.20-45     
##  [79] labeling_0.4.2       rstantools_2.1.1     htmlwidgets_1.5.4   
##  [82] cowplot_1.1.1        bit_4.0.4            tidyselect_1.1.2    
##  [85] processx_3.5.2       plyr_1.8.6           magrittr_2.0.2      
##  [88] R6_2.5.1             generics_0.1.2       DBI_1.1.2           
##  [91] pillar_1.7.0         withr_2.5.0          xts_0.12.1          
##  [94] abind_1.4-5          crayon_1.5.0         car_3.0-12          
##  [97] arrayhelpers_1.1-0   utf8_1.2.2           tzdb_0.2.0          
## [100] rmarkdown_2.13       grid_4.1.3           callr_3.7.0         
## [103] threejs_0.3.3        digest_0.6.29        xtable_1.8-4        
## [106] httpuv_1.6.5         RcppParallel_5.1.5   stats4_4.1.3        
## [109] munsell_0.5.0        bslib_0.3.1          shinyjs_2.1.0